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Abstract. We demonstrate how structured decompositions of unitary operators can 
be employed to derive control schemes for finite-level quantum systems that requirc 
only sequences of simple control pulses such as square wave pulses with finite rise and 
decay times or Gaussian wavepackets. To illustrate the technique, it is applied to 
find control schemes to achieve population transfers for pure-state systems, complete 
inversions of the ensemble populations for mixed-state systems, create arbitrary 
superposition states and optimize the ensemble average of dynamic observables. 
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1. Introduction 

The ability to control quantum-mechanical systems is an essential prerequisite for many 
novel applications triat require the manipulation of atòmic and molecular quantum states 
pj. Among the important applications of current interest are quantum state engineering 
0, control of chemical reactions @, U, ü, 01 , con t r °l of molecular motion ||, selective 
vibrational excitation of molècules || , control of rotational coherence in linear molècules 
||, photo-dissociation ||10|| , làser cooling of infernal molecular degrees of freedom |ÏI, Ï2]|, 
and quantum computation (JT3], [L4], [Ï5|, [Tj], p~7|] . 

Due to the wide range of applications, the immediate aims of quantum control may 
vary. However, the control objective can usually be classified as one of the following: 

(i) To steer the system from its initial state to a target state with desired properties, 

(ii) To maximize the expectation value or ensemble average of a selected observable, 

(iii) To achieve a certain evolution of the system. 

Despite the apparent dissimilarity, these control objectives are closely related. Indeed, 
(D is a special case of @ in which the observable is the projector onto the subspace 
spanned by the target state. (|n|) is a special case of where we attempt to find 
an evolution operator that maximizes the expectation value of the selected observable 
either at a specific target time or at some time in the future. Hence, one of the central 
problems of quantum control is to achieve a desired evolution of the system by applying 
external control fields, and the primary challenge is to find control pulses (or sequences 
of such pulses) that are feasible from a practical point of view and effectively achieve 
the control objective. 

Many control strategies for quantum systems have been proposed. Selective 
excitation of energy eigenstates, for instance, can be achieved using light-induced 
potentials and adiabatic passage techniques [|Ï8|, [L9|, |20|, |2Ï| , which have the advantage 



of being relatively insensitive to perturbations of the control fields and Doppler shifts 
arising from atòmic or molecular motion |23|. Efficient numerical algorithms 
based on optimal control techniques have been developed to address problems such as 



optimization of observables for pure-state [£4], g5], (26| and mixed-state quantum systems 



27, 28]. Quantum feedback control using weak measurements or continuous state 



estimation has been applied to quantum state control problems |29|, [30], [31], [32], [34 
Learning control based on genètic or evolutionary algorithms |35|, [36], |37], [ïï], ^ 
has been a useful tool for quantum control, especially for complex problems for 
which accurate models are not available and in experimental settings |f4T| , [42]. Other 
approaches based on local control techniques or a hydrodynamical formulation [fPfl 
have been suggested as well, and this list is not exhaustive. 

In this paper we pursue an alternative, constructive approach to address the 
problem of control of non-dissipative quantum systems. Note that although real atòmic 
or molecular systems are subject to dissipative processes due to the finite lifetimes of 
the excited states, etc, we can treat these systems as non-dissipative if we ensure that 
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the time needed to complete the control process is significantly less than the relaxation 
times. The technique we develop is based on explicit generation of unitary operators 
using Lie group decompositions. Similar techniques have been applied to the problem 
of controlling two-level systems [Ï6[] and especially partides with spin [pET| , [Í8|j . Here 
we employ decompositions of the type discussed in J49| to derive constructive control 
schemes for iV-level systems. We use the rotating wave approximation (RWA) and 
require that each allowed transition is selectively addressable, for example by applying 
a field of the appropriate frequency, or by appropriate selection rules depending on 
the field polarization. This means we must be able to ensure that each control pulse 
drives a single transition only, and that its effect on all other transitions is negligible. 
These assumptions limit the applicability of this approach to systems for which selective 
excitation of individual transitions is feasible such as atòmic or molecular systems with 
well-separated transition freqüències or partides in anharmonic potentials. Certain 
other factors such as Doppler shifts and inhomogeneous or homogeneous broadening 
must also be taken into account, and may require special consideration in specific 
circumstances. 

However, for systems that satisfy the necessary conditions, the proposed technique 
has some very attractive features. It is constructive and can be used to solve a variety 
of control problems ranging from common problems with well-known solutions such as 
population transfer between energy eigenstates to novel problems such as preparation 
of arbitrary superposition states or optimization of observables for iV-level systems. 
Moreover, although the control schemes derived using this technique depend on the 
effective areas, and to a lesser extent, phases of the control pulses, the pulse shapes are 
flexible, which implies that the control objective can be achieved using control pulses 
that are convenient from a practical point of view such as square wave pulses with 
finite rise and decay times (SWP) or Gaussian wavepackets (GWP). SWP are a realistic 
approximation of bang-bang controls, which play an important role in control theory 
and have been shown to be crucial for time-optimal control Since both SWP and 
GWP can in principle be derived from continuous-wave (CW) làsers using Pockel celis or 
other intensity modulating devices, this also opens the possibility for control of certain 
quantum systems using CW làsers, rather than more complex pulsed làser systems and 
pulse-shaping techniques. 



2. Mathematical and physical framework 

We consider a non-dissipative quantum system with a discrete, finite energy spectrum 
such as a genèric iV-level atom, molecule or particle in an (anharmonic) potential. The 
free evolution of the system is governed by the Schrodinger equation and determined by 
its infernal Hamiltonian H , whose spectral representation is 

N 

H = Y,E n \n)(n\, (1) 

n=l 
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where E n are the energy levels and \n) the corresponding energy eigenstates of the 
system, which satisfy the stationary Schrodinger equation 

H \n) = E n \n), 1 < n < N. (2) 

Although this assumption is not required, we shall assume for simplicity that the 
energy levels E n are ordered in an increasing sequence, E 1 < E 2 < ■ ■ ■ < E N , where 
iV < oo is the dimension of the Hilbert space of the system, and that the eigenstates 
{\n) : n — 1, . . . , N} form a complete orthonormal set. 

The application of external control fields perturbs the system and gives rise to a 
new Hamiltonian H = H + Hj, where Hj is an interaction term. If we apply a field 

f m (t) = 2A m (t) cos(cü m í + <j> m ) = A m (t) [é^<M + e-K^+M] ( 3 ) 

which is resonant with the frequency oj m corresponding to the transition \m) — > \m+l), 
and the pulse envelope 2A m (t) is slowly varying with respect to the frequency u m , then 
the rotating wave approximation (RWA) leads to the following interaction term 

HmUm) = A m {t)d m [e i{ujMm) \m)(m + 1| + e ^ {u)mt+<t}m) \m + l)(m|] (4) 

provided that (a) there are no other transitions with the same frequency u) m and (b) 
off-resonant effects are negligible. Note that the latter assumption is generally vàlid only 
if the Rabi frequency Q m of the driven transition is considerably less than the minimum 
detuning from off-resonant transitions Au m i n , i.e., 

max[fi m (í)] = max[2 A m (t)d m /h] < Auj min , (5) 

where d m is the dipole moment of the transition \m) — * \m + 1). 

The evolution of the controlled system is determined by the operator Ü(t), which 
satisfies the Schrodinger equation 

d ~ í ~ M . 1 ~ 

ïïjU(t) = ÍH + £ H m [f m (t)ÚU(t) (6) 

and the initial condition U (0) = /, where / is the identity operator. 



3. Constructive control using Lie group decompositions 

Our aim is to achieve a certain evolution of the system by applying a sequence of simple 
control pulses. Concretely, we seek to dynamically realize a desired unitary operator 
U (t) at a certain target time t = T. In some cases, we may not wish to specify a target 
time in advance, in which case we attempt to achieve the control objective at some later 
time T > 0. 

To solve the problem of finding the right sequence of control pulses, we apply the 
interaction picture decomposition of the time-evolution operator U(t), 

Ü(t) = É>o(í)#j(í), (7) 
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where Ü (t) is the time-evolution operator of the unperturbed system 

N 

Ü (t) = exp (-iH t/h) = £ e- iE ^ h \n)(n\ (8) 



n=l 



and Üi(t) comprises the interaction with the control fields. To obtain a dynamical law 
for the interaction operator Üi(t), we note that inserting 

ïïjU(t) = H Q Üo(t)Üi(t) + ihÜ (t)jÜj(t) 

M 

HÜ(t) = HoÜoWrit) + iU/m(í)]C>o(í)£j(í) 



m=l 



into the Schrodinger equation (|j) gives 

d . . f M a 1 . . 

= U (ty | ^ F m [/ m (í)]| U (t)Ui(t). 

Applying (|) and the rotating wave approximation Hamiltonian 
leads after some simplification (see [Appendix A| ) to 



d_ 

dt 



Al 



Ui(t) = A m (t)d m /h (x m sm(f) m - y m cos(j) m ) Ui{t) 



m=l 



where we set è m ^ n = \m)(n\ and define 



(9) 

to this equation 
(10) 

(11) 



Hence, if we apply a control pulse f k (t) = 2A k (t) cos(a> m í + 4> k ) which is resonant 
with the transition frequency uj m for a time period tk-i < t < t k (and no other fields 
are applied during this time period) then we have 

Üi(t) = V k (t)Üi(tk-i), (12) 
where the operator V k (t) is 

'd rt 



V k (t) = exp 



h 



A k {t') dt' (x m sin 4> k - y m cos 4> k ) 



(13) 



Thus, if we partition the time interval [0, T] into K subintervals [t k -i, t k ] such that to = 
and tx = T, and apply a sequence of non-overlapping control pulses, each resonant with 
one of the transition freqüències ui m = u a ^, then 



Ü(T) = ?7o( W(T) = e- iA ° T ' h V K V K _ x ■ ■ ■ V u 
where the factors V k are 



\4 = exp 



A k (t) dt [x a{k) sin <f) k - y a{k) cos 4> k 



2A k (t) is the envelope of the kth pulse and a is a mapping from the index set {1, 



(14) 



(15) 
,K} 



to the control index set {1, . . . , M} that determines which of the control fields is active 
for t E [ífc_i,í fc ]. 

It has been shown [4í| that any unitary operator U can be decomposed into a 
product of operators of the type V k and a phase factor e lT = det U, i.e., there exists a 
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positive real number T, real numbers C k and (j) k for 1 < k < K , and a mapping cr from 
the index set {1, . . . , K} to the control-sources index set {1, . . . , M} such that 

Ü = ^ t V k Vk-i···V 1 , (16) 

where the factors are 

V k = exp C k (x a ( k ) sin <p k - y a{k) cos 4> k ) . (17) 

This decomposition of the target operator into a product of generators of the dynamical 
Lie group determines the sequence in which the fields are to be turned on and off. A 
general algorithm to determine the Lie group decomposition for an arbitrary operator 
U is described in [Appendix B . 

Note that in many cases the target operator U is unique only up to phase factors, 
i.e., two unitary operators U\ and Ü% in U(N) are equivalent if there exist vàlues 
n G [0, 2vr] for 1 < n < N such that 

^2 = ^1 (íl^\n){n\^ (18) 

where \n) are the energy eigenstates. For instance, if the initial state of the system is 
an arbitrary ensemble of energy eigenstates 

N 

Po = E w n\n)(n\, (19) 

n=l 

where w n is the initial population of state |n) satisfying < w n < 1 and Y,n=i w n — 1> 
then we have 

t _ fï i \nr,\JB n „„ „-iB n l„ \ \ rVt _ TT ^_7Tt 



^Pof/ 2 T = ^1 fe |n)e ie "«;„e- iy "(n| j C/J = C/^o^í 

i.e., the phase factors e l9n cancel. Thus, if the initial state of the system is an ensemble 
of energy eigenstates, which of course includes trivial ensembles such as pure energy 
eigenstates, then we only need to find a Lie group decomposition of the target operator 
U modulo phase factors, i.e., it suffices to find matrices V k such that 



- / N 

We< 

\n=l 



e w "\n)(n\j = V K V K -i • ■ • V x . (20) 

Note that decomposition modulo phase factors, when sufficient, is more efficient since it 
requires in general up to 2{N — 1) fewer steps than the general decomposition algorithm. 
See |Appendix B| for details. 



4. Choice of pulse envelopes and pulse lengths 

Comparing equations flHp and flT?D shows that 

^ f k A k (t)dt = C k V/c, (21) 

i.e., the effective pulse area of the kth pulse is 2C k where C k is the constant in 
decomposition (|Ï6|) . However, the decomposition does not fix the pulse shapes, i.e., 
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(a) Square-wave pulse (b) Gaussian pulse 




Figure 1. Square wave pulse of length Atk with rise and decay time tq and amplitude 
2Ak (a) and Gaussian wavepacket with q k = and peak amplitude 2Ak (b). 

we can choose pulse shapes that are convenient from a practical point of view such as 
square wave pulses with finite rise and decay times (SWP) and Gaussian wavepackets 
(GWP), which can easily be produced in the laboratory. For instance, in the optical 
regime both SWP and GWP can be created using a combination of continuous-wave 
làsers and Pockel celis or other intensity modulating devices. Moreover, GWP are 
naturally derived from most pulsed làser systems. 



4-1- Square wave pulses 

The pulse area of an ideal square wave pulse of amplitude 2Ak and length Atk is 2AkAtk. 
In order to accurately determine the pulse area of a realistic square wave pulse, however, 
we must take into account the finite rise and decay time r of the pulse. We can model 
the pulse envelopes of realistic SWP [see figure [ï] (a)] mathematically using 

2A k (t) = A k {2 + erf [4(í - r /2)/r ] - erf [4(í - At + r /2)/r ]} (22) 

where erf(x) is the error function 



2 f X 2 

erf(x) = —= / e _í dt. 

V 7T JO 



Although this envelope function may appear complicated, it can easily be checked that 
the area bounded by this function and tj.-i < t < equals the area of a rectangle of 
width Atk — To and height 2A k . Thus, the pulse area f At 2Ak(t) dt of a realistic square 
wave pulse is 2Ak(Atk — To), and equation shows that the amplitude of the pulse 
is determined by 

Ak = aT— x T- x Ck = ïa+ hC \j ' (23) 

where d a ^k) is the dipole moment of the driven transition. 
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To ensure selective excitation, the contribution of Fourier components with Au; > 
Auj min must be negligeable. Noting that the Fourier transform of an ideal SWP (r = 0) 
of length At k and amplitude 2A k is 

/2~sin(±A£ fe Au;) . 

where Au is the detuning from the pulse frequency u m , shows that -F(O) = ^J~^A k At k 
and 

F(Au) _ sin(|Aí fc Au;) 
F(0) ~ ±At k Au ' 

i.e., <C 1 if Atk Au ^ 1- Thus, contributions from Fourier components with 

Au; > Au m i n will be negligible if Aí^ Au;~- n . 

Furthermore, noting that C k < f , the peak Rabi frequency for a square wave pulse 
of length Atk with rise and decay time r is 

2A k (t)d a(k) /h] = -rp^- < ttt^— • (25) 

Hence, the Rabi frequency and the amplitude of the pulse can be adjusted by changing 
the pulse length At k , which allows us to ensure that fl5|) is satisfied, and enforce 
laboratory constraints on the strengths of the control fields. 

We can also give an estimate of the time required to implement arbitrary unitary 
operators given certain bounds on the field strength. If the maximum strength of the 
field produced by the mth làser is A m ^ max , i.e, f m (t) = 2A m (t) cos(u m t + m ) < A m ^ max 
then the time required to perform a rotation by C k on the transition \m) — > \m + 1) 
using a SWP with rise and decay time t is 

Atí = a-^t +To - a — r + r °- (26) 

■ rí m,max u "m ■ rí m,max u 'm 



max 

ífe-i<í<ífc 



[Appendix B| shows that any unitary operator U can be generated up to equivalence (fL8|) 
by performing at most N — m rotations by C < | on each transition \m) —>■ \m + l) for 
m = 1,2, . . . , N — 1. Hence, any unitary operator can be implemented up to equivalence 
using SWP of amplitude A m ^ max in at most time T, where 

N-l N-l / * \ 

T = £ max(Atr P )(iV - m) = £ 1 1~ +r )(N-m). (27) 

m=l m=l \ /i m,m,ax a m ) 

Since two additional rotations on each transition are required to generate U exactly, the 
latter can be accomplished in time T' > J2m=i ma,x(At m WP )(N — m + 2). 

4-2. Gaussian wavepackets 

To model a Gaussian wavepacket [see figure (b)] of peak amplitude 2A k centered at 
t* k = tk-i + \At k) we choose the pulse envelope 

2A k {t) = 2A k exp \-q 2 k (t - At k /2 - t^) 2 ] . (28) 
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The constant qk determines the width of the wavepacket. The pulse area of a Gaussian 
wavepacket is y/7r/q k provided that the time interval At}~ is large enough to justify the 
assumption 

't h - 

In the following we choose q k = 4/ Atk, which guarantees that over 99% of the kth pulse 
is contained in the control interval \t k ^i,t k ] since 

rAt fc /2 



/■í/j r -, Z' + OO 

/ exp -q 2 k (t- Aí fc /2-í fc _!) 2 díw / e 



/ e -# dí = ^-erf(g fc Aí fc /2) 
J-At k /2 q k 

and erf(2) = 0.995322. Thus, © shows that the peak amplitude 2A fe of the GWP is 
determined by 

„ Qk h 4hC k 

A k = —j=z x x C k = —=r- — . (29) 

V 71 " «<r(fc) V n ^kd a (k) 

Again, to ensure selective excitation, the contribution of Fourier components with 
Au> > Au m i n must be negligeable. Noting that the Fourier transform of a Gaussian 
wavepacket with q k = 4/ Atk and amplitude 2A k is 



F(Au>) = ^5 fc exp 



Au; 



F(0) 



4ç, 
ils 

exp(-Au; 2 Aíj;/16) 



At fc A fc 
2^ 



exp(-Aüj 2 A^/16) (30) 



where Au; is the detuning from the pulse frequency u m , shows that 
F(Alu) 



i.e., F j^ty <C 1 if At k Auü ^> 4. Thus, contributions from Fourier components with 
Au; > Au min will be negligible if At k ^> 4Au»~* rr 

Furthermore, noting that C*. < |, the peak Rabi frequency for a Gaussian pulse of 
length At k with = 4/ Atk is 

„_S- t W»] = ^ (31) 

Hence, the Rabi frequency can again be adjusted by changing the pulse length Atk, 
which allows us to ensure that (HD is satisfied and enforce laboratory constraints on the 
strengths of the control fields. 

Again, we can give an estimate of the time required to implement arbitrary unitary 
operators given certain bounds on the field strength. If the maximum strength of the 
field produced by the mth làser is A m ^ max , i.e, f m (t) = 2A m (t) cos(u; m t + 4> m ) < A m , max 
then the time required to perform a rotation by C k on the transition \m) — > \m + 1) 
using GWP with q k = 4/ At k is 

A ,GWP _ 8C! kh < 40Fft /on^ 

^ r m - r- A J - A d ' ( ' 

V ■ rí m,max u "m - rí m,max u 'm 



Since any unitary operator Ü can be generated up to equivalence QT8D by performing at 
most N—m rotations by C k < | on each transition \m) — > |m+l) for m = 1, 2, . . . , N—l, 
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the time required to implement U up to equivalence using GWP of (peak) amplitude 

^m,max ÍS at HlOSt 

T = Emax(At^)(iV - m) = £ f /^. ) - m). (33) 

m=l m=l V^m.maiUm/ 

Since two additional rotations on each transition are required to generate U exactly, the 
latter can be accomplished in time T' > J2m=i niax(Aí^ H/p )(A^ — m + 2). 



5. Physical systems used for illustration 

In the following sections we shall apply these results to various control problems. For 
numerical illustrations of our control schemes, we shall consider 

(i) a four-level model of the electrònic states of Rubidium ( 87 Rb) 

(ii) a four-level Morse oscillator model of the vibrational modes of hydrogen fluoride. 

For Rubidium ( 87 Rb) we consider four electrònic states, which we label as follows: 
|1> = \5S 1/2 ), |2> = |5P 3/2 ), |3) = \AD 1/2 ) and |4) = |6P 3 /2>, where |1) is the ground 
state. Figure || (a) shows the coupling diagram with transition freqüències and dipole 
moments. 



For hydrogen fluoride (HF) we use the Morse oscillator model given in j51| . The 
energy levels corresponding to the vibrational states \n) are 

E n = hüü (n- i) [l - §(n- |) 

where c^o = 0.78 x 10 15 Hz and B = 0.0419. The freqüències for transitions between 
adjacent energy levels are uo n = huo(l — Bn) and the corresponding transition dipole 
moments are d n = poy/n with po = 3.24 x 10~ 31 C m, which leads to the vàlues shown 
in figure (b). Although there are 24 bound vibrational states for this model, we only 
consider the four lowest vibrational modes n = 1, 2, 3, 4, where |1) is the ground state. 

Since we have made several approximations in developing our control approach 
using Lie group decompositions, we must ensure that the assumptions we made are 
vàlid for the systems we consider: 

(i) No two transitions have the same transition frequency.||] 

(ii) Dissipative effects are negligible. 

(iii) The effect of the pulse on off-resonant transitions is negligible. 

Note that both models satisfy hypothesis (p. Furthermore, the main source of 
dissipation for both systems is spontaneous emission. Thus, dissipative effects will be 
negligible provided that the control pulses are much shorter than the lifetimes of the 
excited states. Since the lifetimes of the excited electrònic states for 87 Rb are 28, 90 
and 107 ns, respectively, hypothesis (juD will be satisfied for control pulses in the sub- 
nanosecond regime. Similarly for HF. 

Hypothesis ([ui]) will be satisfied provided that: 

| Assumption (^) can be relaxed if we can distinguish transitions with the same transition frequency 
by other means, e.g., by using fields with different polarizations. 
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(a) 87 Rb 
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(b) HF 
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Figure 2. Transition diagram for Rubidium (a) and hydrogen fluoride (b). For 87 Rb 
the constant u>o = 10 15 Hz and the elèctric dipole moment unit is po = 4.89 x 10~ 29 C 
m. For HF the elèctric dipole moment unit is po = 3.24 x 10 -31 C m. 



(a) The Fourier spectrum of the pulse does not overlap with other transition freqüències, 
i.e., the frequency dispersion of the pulse is less than the minimum detuning from 
off-resonant transitions. 

(b) Equation @ holds, i.e., the Rabi frequency of each driven transition is much smaller 
than the minimum detuning from off-resonant transitions. 

Since the minimum detuning from off-resonant transitions is Au m i n ~ 4 x 10 14 Hz for 
87 Rb and Au min « 3.27 x 10 13 Hz for HF, the pulse length At k should be at least 10" 12 
and 10" 11 seconds, respectively, to ensure that the frequency dispersion of the pulse is 
sufficiently small. Moreover, inserting the vàlues for Auj m i n as well as (p5"D and (|3l|) , 
respectively, into equation (|) shows again that we must choose the pulse lengths such 
that At k > 10" 14 s for 87 Rb and At k > 10" 13 s for HF to ensure that the second 
condition above is met. In the following, we shall choose Atk — 2 x 10" 10 seconds (200 
ps) for all pulses, which ensures that both hypotheses @ and are met for both 
87 Rb and HF. Moreover, such pulses are also experiment ally realizable. 

Note that the energy levels for 87 Rb are multiply degenerate due to hyperfine and 
other effects. Since the detuning between the F = 1 and F = 2 sublevels of the 
55*1/2 ground state is rather large (6.8 GHz), we may wish to be precise and choose 
|1) = \5S\/2,F = 1), for instance, but we shall generally ignore the hyperfine energy 
level structure here. For the cases we consider in this paper, this is justified since the 
frequency differences between the hyperfine levels (except for the ground state) are on 
the order of several hundred MHz or less, which corresponds to detunings of Au; < 10 8 
Hz, which we cannot resolve with 200 ps pulses for reasons outlined above. 
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6. Population transfer |1) — > \N) for a iV-level system 

We shall first apply the decomposition technique described above to the rather 
elementary control problem of population transfer between energy eigenstates to better 
illustrate the technique. Concretely, we consider the problem of transferring the 
population of the ground state |1) to the excited state \N) by applying a sequence 
of control pulses, each resonant with one of the transitions freqüències u> m . It can easily 
be verified that any evolution operator U of the form 



U 






A N -! 


e w 






(34) 



where A^-i is an arbitrary unitary (N — 1) X (N — 1) matrix, e ld is an arbitrary phase 
factor and is a vector whose N — 1 elements are 0, achieves the control objective since 







e i9 



A 



N-l 









and thus the population of state \N) is equal to V e~ ieN e l6N = 1 after application of Ü. 
Next, we observe that setting 

Ü = Ü (T)Üj, Ü I = V N . 1 V N . 2 ···V 1 , (35) 

where the factors are 



V m = exp 



- (x m sin0 m - y m cos« 



— ÍÍp^ m P , 1 4- p-'^ m P n U P 



(36) 



for 1 < m < N — 1, always leads to a U of the form (|34|) , independent of the initial 
pulse phases <p m . 

The factorization (|35f) corresponds to a sequence of N — 1 control pulses in which 
the mth pulse is resonant with the frequency u m of the transition \m) — > \m + l) and has 
effective pulse area tt. Thus, the solution obtained using the decomposition technique 
is an intuitive sequence of 7r-pulses designed to transfer the population step by step to 
the target level. 

The results of illustrative computations for the four-level 87 Rb system introduced 
above are shown in figure [3|. The top graphs show the pulse sequence for square wave 
pulses (a) and Gaussian control pulses (b). The corresponding evolution of the energy- 
level populations shows that the populations of the intermediate levels increase and 
decrease intermittently as expected, while the population of target level |4) reaches one 
at the final time. The bottom graph shows that the energy of the system increases 
monotonically from its kinematical minimum value at t = to its maximum value at 
the final time as predicted. The bàsic response of the system is the same for square 
wave pulses and Gaussian pulses. However, the energy increases more uniformly for 
square wave pulses, while Gaussian pulses tend to result in short, steep increases with 
long intermittent plateau regions. Square wave pulses may therefore be a better choice 
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(a) Square wave pulses (b) Gaussian pulses 




Time(10" 10 s) Time(10" 10 s) 



Figure 3. Population transfcr from the ground state |1) = |55i/2) to the excited 
state |4) = I6P3/2) for 87 Rb using (a) three 200 ps square wave pulses with rise and 
decay time tq — 20 ps and (b) three 200 ps Gaussian pulses with q = 2 x 10 10 
Hz. The top graphs show the pulse envelopes Ak{t). The effective pulse area 
EPA = f At 2A m (t)d m dt of all pulses is ir. The labels 'Field rrí indicate that the 
corresponding pulses are resonant with the frequency u) m of the transition \m) — ► 
|m+l). 

if one wishes to minimize the time the system spends in intermediate states with short 
lifetimes. Gaussian wavepackets, on the other hand, have the advantage of minimal 
frequency dispersion and are thus less likely to induce unwanted off-resonant effects. 

As regards the field strengths, note that for 200 ps pulses up to 380 kV/m are 
required for SWP, and up to 780 kV/m for Gaussian pulses, which corresponds to 
(peak) intensities I = e cE 2 of up to 40 kW/cm 2 (SWP) and 160 kW/cm 2 (GWP), 
respectively. Achieving these intensities experimentally with CW làsers is feasible using 
a combination of sufnciently powerful làsers and beam focusing techniques. Since pulsed 
làser systems with 1 mJ output for picosecond pulses are common, intensities of up to 
10 7 W/cm 2 should be easy to achieve for these systems. 
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Note that we chose pulses of fixed length 200 ps and allowed the pulse amplitudes 
to vary. Had we instead fixed the strength of the fields to be 2Ak = 10 5 V/m, say, 
then the length Atk of the control pulses according to ( pü|) would have been 124.2, 
132.7 and 697.1 ps, respectively, for SWP with r = 20 ps. For Gaussian pulses with 
qk = 4/ Atk, the pulse length according to ( |32"D would have been 235.1, 254.7 and 1528.2 
ps, respectively. Thus, instead of 600 ps in both cases, the time required to achieve the 
control objective would have been 954 ps for SWP and 2018 ps for GWP. 



7. Inversion of ensemble populations for a mixed-state system 



Sequences of 7r-pulses similar to the ones derived in the previous section have played 
an important role in the theory of atòmic excitation [|52| and have been applied to the 
problem of vibrational excitation of molècules in both theory |k| and experiment [54 



The decomposition technique is an important tool since it allows us to generalize the 
intuitive control schemes for population transfer between energy eigenstates to obtain 
similar schemes for a variety of more complicated problems, as we shall demonstrate 
now. 

The first example we consider is a generalization of the population transfer problem 
to mixed-state systems. The objective is to achieve a complete inversion of the ensemble 
populations given an arbitrary initial state of the form (|Ï9"|) . This control operation can 
be regarded as an ensemble-NOT gate for mixed-state systems, not to be confused with 



other NOT-gates such as the U-NOT gate 
populations requires an evolution operator 



Complete inversion of the ensemble 



/ 



U 



V 










,\6 N 



















,i0l 








(37) 



where the e n are arbitrary phase factors. Assuming as before that each transition 
between adjacent energy levels can be individually addressed, the generators of the 
dynamical Lie àlgebra are again of the form ([15]) and the target operator ( pTD can be 
written product of these generators 



ü=ü (t) n 

e=N-i 



UVrr 



.m=l 



(3? 



where the factors V m are as defined in (p6|) . The decomposition ( p8|) corresponds to 
a sequence of K — N(N — l)/2 pulses in which the fcth pulse is resonant with the 
transition \cr(k)) — > \o~(k) + 1) and has effective pulse area ir, where 

a([l,...,K}) = [l,2,···,iV-l;l,2···iV-2;l,2,···,iV-3;···;l,2;l]. 

This pulse scheme does not depend on the vàlues of the initial populations, i.e., a 
complete inversion of the ensemble populations is achieved for any initial ensemble. 
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(a) Square wave pulses (b) Gaussian pulses 




Time(10" 10 s) Time (10 -10 s) 

Figure 4. Inversion of the ensemble populations for the vibrational modes of HF 
using (a) six square wave control pulses with rise and decay time To = 20 ps, and (b) 
six Gaussian pulses with q = 2 x 10 10 Hz. The top graphs show the pulse envelopes 
Ak(t). The effective pulse area EPA = J At 2A m (t)d m dt of all pulses is ir. The labels 
'Field rrí indicate that the corresponding pulses are resonant with the frequency uj m 
of the transition |m) — > |m + 1). 

Moreover, if the initial populations are mutually distinct, i.e., w n ^ w m for n ^ m, then 
the decomposition is optimal in the sense that a complete inversion of the ensemble 
populations cannot be achieved with fewer than K control pulses. 

To illustrate the control scheme, let us apply it to the four-level Morse oscillator 
model for the vibrational modes of HF discussed above. For the purpose of the computer 
simulations, we randomly choose the initial populations to be w\ = 0.4, u> 2 = 0.3, 
w 3 = 0.2 and w 4 = 0.1, but recali that any initial ensemble would do, i.e., we 
could have chosen a thermal ensemble given by a Boltzmann distribution or another 
ensemble instead. Our goal is to create an ensemble where the populations of the 
energy eigenstates are reversed, i.e., where |1) has population W4, |2) has population u> 3 , 
|3) has population w 2 , and |4) has population w\. 
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Figure fl shows the results of control simulations using square wave and Gaussian 
control pulses, respectively. Note that each pulse in the control sequence interchanges 
the populations of two adjacent energy levels until a complete inversion of the 
populations is achieved. For our four-level system the effect of the controls on the 
populations can be summarized as follows 

fi h h fi h fi 
\1) Wi ,w 2 w 2 w 2 ,W 3 w 3 w 4 

\2) w 2 wi w 3 w 3 w 2 Wi 

|3) w 3 w s Wi W4 W4 ^w 2 Wi 

|4) w 4 W4 W4 w± wi wi wi 

where f m , m = 1,2, 3, refers to a control pulse of frequency u) m with effective pulse area 
n. The first pulse interchanges the populations of levels |1) and |2), the second pulse flips 
the populations of levels |2) and |3), the third pulse switches the populations of levels |3) 
and |4), etc. Since the populations of our initial ensemble satisfy Wi < w 2 < w 3 < w 4 , 
the energy of the system assumes its kinematical minimum at t — and increases 
monotonically to its kinematical maximum value at the final time. Again, the gradient 
of approach is more uniform for square wave pulses. 

As regards the field strengths, note that for 200 ps pulses up to 5.7 MV/m are 
required for SWP, and up to 12 MV/m for Gaussian pulses, which corresponds to 
(peak) intensities / = e cE 2 of up to 8.5 MW/cm 2 (SWP) and 24 MW/cm 2 (GWP), 
respectively. Achieving these intensities experiment ally should be no problem for pulsed 
làser systems. For CW làsers, it would be challenging at the moment, but it should 
still be feasible using a combination of powerful làsers and beam focusing techniques. 
Moreover, such problems should disappear with improvements in làser technology in the 
future. 

Had we instead of fixing the pulse length at 200 ps, fixed the strength of the fields 
to be 2Ak = 5 x 10 6 V/m, say, then the length At^ of the control pulses according to (^6|) 
would have been 224.5, 164.6, 138.1, 224.5, 164.6 and 224.5 ps, respectively, for SWP 
with r = 20 ps. For GWP with q k = 4/Aí fc the pulse lengths according to fl32|) would 
have been 461.3, 326.2, 266.3, 461.3, 326.3 and 461.3 ps, respectively. Thus, instead of 
1.2 ns in both cases, the time required to achieve the control objective would have been 
1.14 ns for SWP, and 2.3 ns for GWP. 

Note that the problem of population transfer for a system initially in state |1) is a 

special case of the problem of population inversion for a trivial ensemble with w± — 1 

and w 2 = w 3 = w 4 = 0. It can easily been seen that pulses four, five and six in the pulse 

sequence above do no harm but have no effect for this initial ensemble 

fi fï Í3 fi f'2 fi 
|1) 1 

|2) ^1 ^0 ^0 

|3) 'M ^0 1 

|4) ^1 1 1 1 
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and can therefore be omitted. Thus, the general six pulse sequence simplifies in this case 
to the three pulse sequence in the previous section. This can also be inferred directly 
from the decomposition (]3l|) of the target operator. For a four level system ([38]) becomes 
U = U (T)ViV2ViV 3 V2Vx with V m as in (|36|). Thus, after applying the pulse sequence 
the state of the system is 

p = ^o(T)yií> 2 í>iy 3 ^t>iPo[^o(r)yiWiHV r 2^i] t 

If p = |1)(1| then V 3 V 2 ViPoV?Vlvj = |4)(4|. Since Ü (T)ViV 2 Vi commutes with this 
operator, the remaining factors cancel in the decomposition 

= f>o(T)V 1 V 2 V 1 |4)(4|V 1 t V r 2 t V 1 t f> (T) t 
= |4)(4|f> (T)f 1 y 2 V r 1 V r 1 t y 2 t Vl t ?7o(T) t = |4)(4| 
and hence do not change the state of the system. 



8. Creation of arbitrary superposition states 



In this section we consider the problem of creating arbitrary superposition states from 
an initial energy eigenstate. Control schemes to create such superposition states may be 
useful in controlling quantum interference in multi-state systems, and can be considered 
a generalization of the ~ pulses used routinely in free induction-decay experiments [56]. 

Concretely, assume that the system is initially in the ground state |1). To create 
the superposition state 



N 



!*(*)> 



J2r n é 9 -é Ent/h \ 



n) 



n=l 



N 
n=l 



n{t)) 



(39) 



where the coefficients r n satisfy the normalization condition Y,n=i r n = 1> we need to 
find a unitary operator Üj such that 

N 



Ut\1) = Y j r n é e -\n) 



(40) 



n=l 



and decompose Uj according to the algorithm described in [Appendix B 
To find a unitary operator Üj that satisfies (B4T) we set 



W 



r-2 



\ 



'Af-l 



(41) 



where Ijv-i is the identity matrix in dimension N — 1, and perform Gram-Schmidt 
orthonormalization on the columns of W. This produces a matrix Ü\ which is unitary 
and satisfies £7i|l) = Hn=i r n\ n )- Hence, Üi = QÜi with 6 = J2n=i e l9n \n){n\ satisfies 
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As an example, we consider the problem of creating the superposition state 
I^W) = f£ n =i I^W) f° r a four-level system initially in state |1). As outlined above, 



we set 



W 



( 1/2 \ 

1/2 10 

1/2 10 

V 1/2 1 ) 



(42) 



Ui 



(43) 



and perform Gram-Schmidt ort honor malizat ion on the columns of W, which gives 

/ 1/2 -V3/6 -VE/6 -y/2/2 \ 

1/2 +y/3/2 

1/2 -^3/6 +^3 

V 1/2 -V3/6 -^6/6 +^2 / 

Since = / we have Üj = Ü\ and applying the decomposition algorithm (appendix 
Appendix B|) leads to the factorization Ü\ = V5V4V3V2V1, where the factors are 



Vi 


= exp (+CxXi) , 


C x = 


TT 

3' 




= exp (-C 2 x 2 ) , 


c 2 = 


arctan (\/2) 


% 


= exp (+C3X3) , 


c 3 = 


7T 

4 " 


Va 


= exp (+C4X2) , 


c 4 = 


TT 

2 ' 


% 


= exp (-C5X1) , 


c 5 = 


TT 

2 ' 



(44) 



This decomposition corresponds to the following sequence of five control pulses 





= A^é^ 


-7T/2) 


+ C.C. 


= -2A x {t) singlí) 


/ 2 (í) 


= A 2 (t)e i ^ 2t - 


-tt/2) 


+ C.C. 


= +2A 2 {t) ún{u 2 t) 


/s(f) 


= Az{t)é^ zU 


-tt/2) 


+ C.C. 


= -2A 3 (t) sm(üü 3 t) 


k{t) 


= A 4 (í)e í ^ 2ÍH 


-tt/2) 


+ C.C. 


= -2A 4 (t) ún(uj 2 t) 




= A 5 (t)e i( - Wlt - 


-tt/2) 


+ C.C. 


= +2A 5 (t) singlí) 



with pulse areas |7r, 2 arctan(\/2), ^n, tt and tt, respectively. Note that only five instead 
of six pulses are required since the target operator U\ has two consecutive zeros in the 
last column, which implies that one of the six control pulses has zero amplitude and can 
thus be omitted. 

Figure [| shows the results of a control simulation based on this decomposition of 
U\ for 87 Rb using square wave and Gaussian control pulses, respectively. Note that all 
the populations and the absolute vàlues of all the coherences are 0.25 at the final time 
- exactly as required for the superposition state \^{t)) = | J2n=i whose density 

matrix representat ion is 



1 

-iüji2í 

-ÍOJ14Í 



-ÍW23* 



3 Í1^23Í 



e iwi 4 t \ 
D Íí^24Í 



3 ÍW34Í 



g— ÍW24Í g~ 1^34* 
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(a) Square wave pulses (b) Gaussian pulses 




02468 10 02468 10 



Time (10" 10 s) Time (10 -10 s) 

Figure 5. Creation of the superposition state | \& (t) ) = | J2n=i l^-(O) f° r 87 Rb initially 
in the ground state |1) using (a) fïve 200 ps square wave pulses with rise and decay 
time To = 20 ps, and (b) five 200 ps Gaussian pulses with q = 2 x 10 10 Hz. The top 
graphs show the pulse envelopes Ak(t). The effective pulse area (EPA) of all pulses is 
as shown in the graph. The labels 'Field m' indicate that the corresponding pulses are 
resonant with the frequency uj m of the transition |m) — > |m + 1). 



i.e., \p mn \ = i for all m, n. Note that we have plotted the absolute vàlues of the 
coherences p mri (t) (for m ^ n) only since their phases are rapidly oscillating at 
freqüències u mn = (E n — E m )/h, which are on the order of 10 15 Hz for 87 Rb. The 
pulse intensities are similar to those for population transfer in 87 Rb. Again, we chose 
pulses of fixed length 200 ps. Had we instead fixed the strength of the fields to be 
2Ak = 10 5 V/m, say, then the length At^ of the control pulses according to (26) and 
(j32|) would have been 99.5, 88.6, 358.6, 132.9 and 155.4 ps, respectively, for SWP with 
r = 20 ps, and 156.7, 154.9, 764.1, 254.7 and 305.6 ps, respectively, for GWP with 
q k = 4/Aí fc . 

Unlike decompositions (|35|) and (|38D where the initial phases m of the control 
pulses were arbitrary, the factorization (041) fixes the pulse area and frequency u m as 
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well as the initial phase m of each pulse. In order to determine the significance of the 
the initial pulse phases on the outcome of the control process, we compute the unitary 
operator U 2 = \ A \\ A A 1 . where the factors are 

V\ = exp [Ci (singlí 1 — cos fayi)] 
V 2 = exp [C 2 (sin fax 2 - cos fam)] 

V 3 = exp [C 3 (sin03X3 - cos0 3 ^ 3 )] (45) 

V4 = exp [C 4 (sÍn 04^2 — cos 4 y 2 )] 

V 5 = exp [C 5 (sin 05^1 - cos 5 yi)] 

and the constants are as in (fP]) but the initial phases 0& of the control pulses are 
arbitrary, and apply this operator to the initial state The resulting state 



U, 



( 1 ^ 




/ 


gÍ((^4+</'5-</>l-^2) 


\ 





1 




e i(-ir/2-0 5 ) 







~ 2 




gí( 7r— <t>i—4>^) 




V j 




V 


gi(7r/2— <^l— <f>2-<t>3) 


) 



(46) 



differs from the desired target state only in the phase factors, i.e., the pulse phases do 
not affect the relative amplitudes r n of the superposition state created. Moreover, we 
can use (ffl) to explicitly determine the pulse phases n as a function of the phases 9 n 
of the target state: 

arbitrary 

| - 20i - 0! - 2 - e 3 

0i + 9 1 + 9 2 + 9 3 - 4 (47) 

7T - 0i - 3 

1-9, 



(pi - 

02 = 

03 = 

04 = 

05 = 

Setting 9 n = for n 
04 = vr/2 and 05 = - 



2 ^2- 

= 1,2,3,4 and choosing 0i = n/2 leads to 02 = — tt/2, 3 
•7r/2, which agrees with the phases in decomposition (fHj). 



vr/2, 



9. Optimization of observables 

Finally, we address the problem of maximizing the ensemble average of an observable 



for a system whose initial state is a statistical ensemble of energy eigenstates (19). Let 
us first consider the case of a time-independent observable A. To determine the target 
operator required to maximize the ensemble average (A) of A we observe that (A) is 
bounded above by the kinematical upper bound 



N 

(i) < £«V(„)A n , (48) 

n=l 

where À n are the eigenvalues of A counted with multiplicity and ordered in a non- 
increasing sequence 

Ai > À 2 > ••• > A^v, (49) 



Constructive control of quantum systems usíng factorization of unítary operators 21 



w n are the populations of the energy levels E n of the initial ensemble, and o~ is a 
permutation of {1, . . . , N} such that 

w a (i) > w a{ 2) >■■■> w a{N) . (50) 

Observe that this universal upper bound for the ensemble average of any observable 
A is dynamically attainable since the systems considered in this paper are completely 
controllable [|^, 

Let |\I/ n ) for 1 < n < N denote the normalized eigenstates of A satisfying 
= A n |\I/ n ) and let Üi be a unitary transformation such that 

!*„(„)> = Üiln), l<n<N. (51) 

Given an initial state po of the form (|Ï9|), we have 



Ti(ÀÜ lPo Ül) = Tr\À^w n Ü 1 \n)(n\ÚU 

= Tr ( J2 W n À\ }(* CT (n)| 



J2w n X a ( n ) = ^2w a ( n )X n . (52) 



Hence, if the system is initially in the state (|Ï9|) then U\ is a target operator for which 
the observable A assumes its kinematical maximum, and we can use the decomposition 
algorithm described in |Appendix B| to obtain the required factorization of the operator 

üi = üoiryü,. 

However, if A is an observable whose eigenstates are not energy eigenstates then the 
expectation value or ensemble average of A will usually oscillate rapidly as a result of 
the action of the free evolution operator Uo(t). These oscillations are rarely significant 
for the application at hand and often distracting. In such cases it is advantageous to 
define a dynamic observable 

À{t) = Ü (t)ÀÜ (ty (53) 

and optimize its ensemble average instead. To accomplish this, note that if \^ n ) are 
the eigenstates of A satisfying A|\l> n ) = X n \^f n ) then \^/ n (t)} = U (t)\^ n ) are the 
corresponding eigenstates of A(t) since 

À(t)\Ü n (t)) = Ü Q (t)ÀÚ (tfU (t)\V n ) = Ü (t)\ n \V n ) = An|*„(í)>. 

Hence, if U\ is a unitary operator such that equation (|5TD holds then Uç,(t)Ui is a unitary 
operator that maps the energy eigenstates |n) onto the A(í)-eigenstates \^ n (t)) since 

Üoitp^n) = £/„(*) | = |*a(n)(*)) 

for 1 < n < N. Thus, the evolution operator required to maximize the ensemble 
average of A(t) at time T > is U (T)Ui and the target operator to be decomposed is 
U = Ü (T)iÚ Q (T)Ü 1 = Ü 1 . 
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For instance, suppose we wish to maximize the ensemble average of the transition 
dipole moment operator À(t) = U (t)ÀÜ (ty , where 

2V-1 

À=^d n {\n)(n + l\ + \n + l)(n\) } (54) 

n=l 

for a system initially in state (O) with 

w\ > w 2 > ■ ■ ■ > wn > 0. (55) 

First, we need to find a unitary operator that maps the initial state \n) onto the A- 
eigenstate |\l/ n ) for 1 < n < N. Let U\ be the N x N matrix whose nth column is the 
normalized A-eigenstate \^ n )- Then U\ clearly satisfies U\\n) = \^ n )- Furthermore, XJ\ 
is automatically unitary since the eigenstates |\l/ n ) are orthonormal by hypothesis. 

For N = 4 and d n = poy/n the eigenvalues of the operator A defined in (|54]) are (in 
decreasing order) 



Ai = V3 + \/6, A 2 = V3 - VQ, A 3 = -A 2 , A 4 = -Ai 

and the corresponding eigenstates with respect to the Standard basis \n) are the columns 
of the operator 



Ui 
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2Ai 


2A 2 


2A 2 


2Aj 
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2 
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2 




V2-V3 


V2-V3 


V2+V3 


2Ai 


2X 2 


2X 2 


2Ai 
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2 


2 


2 



(56) 



Applying the decomposition algorithm described in |Appendix B| yields the product 
decomposition ÜiQ = VqVVVVV, where the factors are 

tt/4, 

arctan (v^) , 



vi 


= exp 




d 


v 2 


= exp 


[-C2X2) , 


c 2 


% 


= exp 


'-C3X1) , 


c 3 


Va 


= exp 


^-C 4 x 3 ) , 


c 4 


% 


= exp 


[-C 5 x 2 ) , 


c 5 


V e> 


= exp 


{-C 6 xi) , 


c 5 



3 + 



tt/3, (57) 
arctan 
arccot 

and 6 = diag(l, —1, 1, —1). Note that Ü 2 = Éi© is equivalent to since 9 commutes 
with po as defined in equation (|Ï9D, i.e., 0p o 0^ = po, and thus 

Tr {ÀÜ2P0ÜÍ) = Tr (ÀÜxèpo&ÚÏ) = Tr (itfipWi) . (58) 

This decomposition corresponds to a sequence of six control pulses 



fi(t) 

h{t) 
hit) 
h{t) 
h(t) 

fe(t) 



A 1 (t)e i( - Wlt - W ^ + c.c. 
A 2 {t)é^ 2t - n ^ + c.c. 
Azity^-^W + c.c. 
Aiity^-^W + c.c. 
Asity^^W + c.c. 
Ael^e^-^ + c.c. 



2A 1 (t) singlí) 
2A 2 (í) sin(u; 2 í) 
2A 3 (í) singlí) 
2A 4 (í) sin (u; 3 t) 
2A 5 (í) sin(u; 2 í) 
2A 6 (t) sin(o;ií) 



Constructive control of quantum systems usíng factorization of unítary operators 23 



with effective pulse areas |-, 2C 2 , 2C 3 , 4f , 2C 5 and 2C 6 , respectively. A gain, the 



3 • 



decomposition fixes the frequency and pulse well as the initial phase of each pulse 

and the question thus arises what role the phases play. As we have already seen, the 
target operator U\ is not unique. In fact, equation (|58| ) shows that right multiplication 
of U\ by any unitary matrix that commutes with po produces another unitary operator 
that leads to the same ensemble average of the target observable. Nevertheless, in 
general, the control process is sensitive to the phases m . For instance, one can verify 
that changing the phase 0i of the first pulse from — 7r/2 to 7r/2 in the pulse sequence 
above leads to the following evolution operator 



1 


1 
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2A 2 
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V2+-/3 
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-2Ai 
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(59) 



which maps |3) onto and |4) onto — ^4) but |1) onto |\I/ 2 ) an d |2) onto and 
leads to the ensemble average 

(A) = Wi\ 2 + w 2 \i + W3À3 + u; 4 A 4 (60) 

at the final time, which is strictly less than the kinematical maximum if w± > w 2 . 

Figure H shows the results of control simulations for HF with initial populations 
Wi = 0.4, w 2 = 0.3, W3 = 0.2 and w 4 = 0.1 for square wave and Gaussian control pulses, 
respectively. The pulse intensities are similar to those for population inversion in HF. 
Notice that the observable indeed attains its kinematical upper bound at the final time, 
as desired. Furthermore, the target state for which the observable assumes its upper 
bound is 
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with p l2 = AiA 2 (A 2 + Ai/3)/40 w 0.0658, p u = AiA 2 (A 2 - Ai/3)/40 w -0.0016, 
p 23 = A 2 (A? - 1/^/40 w 0.0904, p u = A 2 (A? + 1/^/40 w 0.1118, which agrees 
with the final vàlues of the populations and coherences in figure ||. Note that we chose 
pulses of fixed length 200 ps. Had we instead fixed the strength of the fields to be 
2Ak = 5 x 10 6 V/m, say, then the length At^ of the control pulses according to (p6|) 
and © would have been 122.2, 97.9, 60.8, 156.3, 82.5 and 72.7 ps, respectively, for 
SWP with r = 20 ps, and 230.6, 198.4, 92.2, 307.5, 141.0 and 118.9 ps, respectively, 
for GWP with q k = 4/Aí fc . 



10. Conclusion 



We have presented several control schemes designed to achieve control objectives ranging 
from population transfer and inversion of ensemble populations to the creation of 
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(a) Square wave pulses 



(b) Gaussian pulses 




4 6 

Time(10" 10 s) 



Time (10" lu s) 



Figure 6. Maximization of the transition dipole moment for HF using (a) six square 
wave pulses with rise and decay time t = 20 ps, and (b) six Gaussian pulses with 
q = 2 x 10 10 Hz (right). The top graphs show the pulse envelopes Ak(t). The vàlues of 
the constants CV- which determine the effective pulse areas (EPA) are given in (j57|) . The 
labels 'Field rrí indicate that the corresponding pulses are resonant with the frequency 
Lu m of the transition \m) — > \m + 1). 
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arbitrary superposition states and the optimization of (dynamic) observables. A key 
feature of these schemes is that they rely only on sequences of simple control pulses 
such as square wave pulses with finite rise and decay times or Gaussian wavepackets 
to achieve the control objective. In the optical regime, for instance, such pulses can 
easily be created in the laboratory using pulsed làser sources, or by modulating the 
amplitude of CW làsers using Pockel celis. No sophisticated pulse shaping technology 
is required. A limitation of the approach is the need to be able to selectively address 
individual transitions, which restricts the application of this technique to systems where 
selection rules and frequency discrimination can be employed to achieve this. However, 
these requirements can be met for certain atòmic or molecular systems, as we have 
demonstrated for Rubidium and hydrogen fluoride. 
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Appendix A. Derivation of equation (|Ï0|) 

Let E n = E n /h and d n = d n /%. Inserting equations (|§p and (Q) into (|Sp leads to 

ü (ty { J2 H m [f m (t)]/n) üoitpjit) 

{m=l J 

^ ] 6 è nn A m {t^d m ( 6 ^ ^ èm,m+í& ^ ^ ^è m _j_x m ^ 6 è n i ^ n 'tj i(t) 

n,m,n' 

= Y^ A mit)d m (e i · Èmí e i(üJraí+ ^ ) e- iÉ '" +lí è mim+ ie iÉm + lí e- i{üJmí+0m) e- iÉmí è m+ i, m ) Ü^t) 

m 

^ ^ A m (t)d m ^6 è m ^ m +\ -\- 6 è m _|_i im ^ t//(í) 

m 

= ^2A m (t)d m [cos 4> m (è m>m+ i + è m+ i jm ) + isin0 m (é m , m+ i - è m+ i, m )] [//(*) 

m 

= J2 A m(t)drn ( Win COS (j) m + X m Í SÍïl m ) £//(í). 

m 

Hence, multiplying both sides by —i gives 

—T7— = J2A m (t)d m (x m sin0 m - y m cos(f) m ) Ü^t). (A.l) 
at m 

Appendix B. Lie group decomposition algorithm 



.düijt) 
1 dt 



To find a decomposition QTB|) for the unitary operator U we define the equivalent operator 
Ü<® G SU(N) by Ü<® = e-' lT l N Ü where e ir = det(£>). Our goal is to reduce Ü<® step 
by step to a diagonal matrix whose diagonal elements are arbitrary phase factors e 10n . 
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Recali that this reduction is always sufficient if the initial state of the system is an 
ensemble of energy eigenstates. 

Let Uif 1 denote the ïth row and jth column entry in the matrix representation of 

. In the first step of the decomposition we seek a matrix 

= exp [—Ci (sin0iXi — cos0i?/i)] , 

which is the identity matrix everywhere except for a 2 x 2 block of the form 

/ cos(Ci) ie^ 1 sin(Ci) \ 
^ ie- 1 * 1 sin(Ci) cos(Ci) ) 

in the top left còrner, such that 
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where c is some complex number. Noting that U[ N 
easily be verified that setting 

C fc = -arccot(-r 2 /ri), <f) k = ir/2 + a x - a 2 

). Next we set Ü w = W^Ü (0) and find W {2) of the form 

= exp [— C 2 (sin0 2 ^2 - cos fom)] 



achieves 
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where c is again some complex number. Repeating this procedure N — 1 times leads 
to a matrix U^" 1 ' whose last column is (0, . . . , 0, e ldN ) T . Since we are not concerned 
about the phase factor e l0N in this paper, we stop here. Note that 

exp (— Cxn-i) x exp [— C(xat_i sin (f> — í/n-i cos </>)] 

with C = 7r/2 and = — tt/2 — 6 n maps (0, e" 9 ^- 1 ) 71 onto (0, 1) T . Hence, a complete 
reduction to the identity matrix would require two additional steps to eliminate e 10N , 
which would result in two additional control pulses. 

Having reduced the last column, we continue with the (N — l)st column in the 
same fashion, noting that at most N — 2 steps will be required to reduce the (N — l)st 
column to (0, ... , 0, e" 9 ^- 1 , 0) T since is unitary. We repeat this procedure until after 



at most K = N(N — l)/2 steps is reduced to a diagonal matrix diag(e 1É>1 , 
and we have 



D i9 N 



WW...WWÚ<® =diag (í 



i0l 



)■ 



(B.7) 



Constructive control of quantum systems using factorization of unitary operators 27 

Finally, setting V k = ÍW^ K+1 - k ^ leads to 

Ü {0) = VkVk-i • ■ ■ V^idiag (e i91 , . . . , e ie -) (B.8) 

and therefore Ü = VrVk-i • • • where 9 = e ir//7V diag (e l8l } . . . , e 10JV ) is a diagonal 
matrix of phase factors. 

Recali that U can always be decomposed such that 8 is the identity matrix. 
However, up to 2(N — 1) additional terms would be required to eliminate the phase 
factors, which would result in additional control pulses. While some applications indeed 
require the elimination of these phase factors, they are often insignificant and the 
additional control pulses would be superfluous. For a more sophisticated decomposition 
algorithm that requires only very few phases the reader is referred to pü . 
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